Global structure of Choptuik's critical solution in scalar field collapse 
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At the threshold of black hole formation in the gravitational collapse of a scalar field a naked 
singularity is formed through a universal critical solution that is discretely self-similar. We study 
the global spacetime structure of this solution. It is spherically symmetric, discretely self-similar, 
regular at the center to the past of the singularity, and regular at the past lightcone of the singularity. 
At the future lightcone of the singularity, which is also a Cauchy horizon, the curvature is finite and 
continuous but not differentiable. To the future of the Cauchy horizon the solution is not unique, 
but depends on a free function (the null data coming out of the naked singularity). There is a unique 
continuation with a regular center (which is self-similar). All other self-similar continuations have 
a central timelike singularity with negative mass. 



I. INTRODUCTION 

In general relativity, a black hole may be formed during 
the evolution from asymptotically flat initial data where 
none is present. Consider a one-parameter family of reg- 
ular asymptotically flat initial data. It is not difficult 
to find such families which form a black hole for some 
range of the parameter (strong data) but disperse for an- 
other range (weak data) . The boundary between the two 
regimes is the black hole threshold. In what is called type 
II critical collapse, the black hole mass can be made ar- 
bitrarily small by adjusting the parameter p of the initial 
data to its critical value p». Near the threshold, the final 
black hole mass M then scales as 



(1) 



where C is a constant. C depends on the family, but the 
transcendental number 7 is universal - it depends on the 
type of matter but not on the family of initial data. 

Type II critical collapse was originally discovered by 
Choptuik in the spherically symmetric massless scalar 
field , but has since been found in many simple matter 
systems in spherical symmetry, and also in axisymmetric 
gravitational waves Q • A review of the field is Q . 

Type II critical phenomena can be described in dy- 
namical systems terms: the phase space of the system 
has (at least) two attracting fixed points, namely black 
holes and dispersion. The boundary between the two 
basins of attraction, the critical surface, contains a crit- 
ical point: it is an attractor within the boundary sur- 
face, and a repeller only in the one remaining direction. 
This means that it must have precisely one unstable lin- 
ear perturbation, with the property that adding a bit of 
that perturbation with one sign leads to collapse, while 
adding it with the opposite sign leads to dispersion. In 
type II critical collapse the critical point is either a dis- 
cretely self-similar (DSS) or a continuously self-similar 
(CSS) spacetime. 



Type II critical collapse is interesting, among other 
things, because the maximum value of the curvature in a 
subcritical evolution, and the maximal value of the cur- 
vature outside the black hole in a supercritical evolution, 
both diverge as 
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as the fine-tuning is improved. (This is basically the same 
result as the black hole mass scaling, and similar results 
hold for any curvature invariant.) From the dynamical 
systems picture it is clear that the end point of type II 
critical collapse in the limit of perfect tuning of p to its 
critical value is not a "zero mass black hole" but the 
critical solution itself. This solution has a naked singu- 
larity. It is therefore interesting to examine the global 
spacetime structure of the critical solution, and in par- 
ticular its Cauchy horizon. Here we do this for the spher- 
ically symmetric massless scalar field, where the critical 
solution is DSS. We focus on this system because CSS 
can be viewed as a limiting case of DSS, and because the 
critical solution in the most interesting system in which 
type II critical phenomena have been found, axisymmet- 
ric pure gravity, is also DSS Q- 

In Section [H] we discuss the global structure of Chop- 
tuik's critical solution kinematically. Section ITTT1 sets out 
the field equations for the real massless scalar field in 
spherical symmetry, in coordinates adapted to our prob- 
lem, and describes the mathematical structure of the so- 
lution at the Cauchy horizon of the singularity. Sections 
II VI IV! and IV1I show the results of our numerical integra- 
tion of the critical solution, and Section fVIII contains our 
conclusions. Some technical details have been removed 
to appendices. 



II. KINEMATIC AL DISCUSSION 

A spacetime is discretely self-similar (DSS) if there is 
a conformal isometry $ of the spacetime such that 
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The value of the dimensionless "logarithmic scale period" 
A is a geometric property of the spacetime, independent 
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of coordinates. It is often useful to work in coordinates 
adapted to the symmetry at hand. A generic self-similar 
and spherically symmetric metric can be written as 

ds 2 = e~ 2T (A dr 2 + 2Bdrdx + C dx 2 + F 2 dfl 2 ) (4) 

where dQ 2 = d9 2 + sin 2 9 dip 2 is the metric on the unit 
2-sphere, and where A, B, C and F are functions of r 
and x. This metric is DSS if and only if they are periodic 
in r with period A. It is continuously self-similar (CSS) 
if these functions are completely independent of r. We 
assume that the signature is (—,+,+,+), and that the 
metric is non-degenerate. This leads to the inequality 
AC — B 2 < 0. We also assume F > for definiteness. 

Any four-dimensional spacetime splits into a product 
of a two-dimensional spacetime (the reduced manifold) 
and a round two-sphere of area Anr 2 . The area radius 
r is a scalar in the reduced manifold. Here the coor- 
dinates on the reduced manifold are r and x, and the 
area radius is given by r — e~ T F. Geodesies in the re- 
duced spacetime are radial (constant 9 and ip geodesies 
in the full spacetime). The Hawking mass m is defined 
by 1 — 2m /r = (Vr) 2 . It is a scalar on both the full 
and the reduced spacetime. From m we define the two 
dimensionless scalars [l = 2m/r and a = (1 — /i) -1 / 2 . 
It is easy to show that a spherical surface where /Lt > 1 
is a closed trapped surface, and one where fi = 1 is an 
apparent horizon. In a DSS spherical spacetime, \i and 
a are periodic in r. 

Radial null geodesies which are invariant under the 
symmetry @ are called self-similarity horizons (SSH). 
They are the key to determining the causal structure. 
All coordinate systems of the form (QJ are related by co- 
ordinate transformations of the form 

x' = ip{r,x), t' = r + ip(r,x), (5) 

where if and tp are periodic in r with period A. We use 
this coordinate freedom to make all lines in the reduced 
manifold where F — into lines of constant x. (These 
can be either regular centers or central singularities.) We 
also make all SSHs into lines of constant x where A = 0. 

In order to discuss the global structure of the Choptuik 
spacetime and its possible continuations we briefly review 
the kinematical results of 4] . In a spherically symmetric 
DSS spacetime, two kinds of singularities can be distin- 
guished. From dimensional analysis it can be seen that 
the Kretschmann scalar scales as e 4r for constant x, and 
therefore the set r = oo is a central (because r scales 
as e~ T ) curvature singularity. We call this the kinemati- 
cal singularity. Geometrically, this singularity is either a 
point or a null line in the reduced spacetime. There are 
two types of self-similarity horizons that in £| we have 
called fans and splashes. The kinematical singularity is 
null if there is at least one splash. Additional central 
singularities can arise where F = for all r. (We call 
these dynamical.) Because r takes all values up to oo 
they are connected to the kinematical singularity. There 
are at most two of them, connected to the kinematical 



singularity at its ends. Topologically, they are lines in 
the reduced manifold. 

All known type II critical solutions in spherical gravi- 
tational collapse can be defined by the properties of self- 
similarity (CSS or DSS), analyticity at the past light- 
cone, and the requirement that they have a single unsta- 
ble perturbation mode. A generic spherically symmetric 
DSS scalar field solution is singular either at the center 
or the past lightcone. Imposing analyticity at both the 
center and the past lightcone defines a nonlinear PDE 
boundary value-problem which admits at most discrete 
solutions 

Only one such solution has been found, and 
it empirically turns out to have only one unstable mode, 
and to agree perfectly with the critical solution found 
previously in collapse simulations by Choptuik [lj. 

The global structure of the Choptuik solution up to the 
future light cone of the kinematical singularity (which is 
a Cauchy horizon) is sketched in Fig. 2] together with 
the x and r-lines of the three coordinate patches that we 
shall use in the numerical calculations. This structure is 
the same as for all other known type II critical solutions 
in spherical symmetry. These solutions have a regular 
center x — x c in the past. As x increases, the x lines 
are at first timelike, so that A < 0. They become null 
at the fan x = x p where ^4 = 0, dA/dx > and B > 
0. As x increases further, they are spacelike, so that 
A > spacelike. Somewhere in the spacelike region B 
changes sign. The s-lines become null again at the second 
fan x — Xf, where A = 0, dA/dx < and B < 0. 
Approaching the kinematical singularity t — oo from the 
range of "angles" x c < x < Xf it is a single point. x p 
and Xf are its past and future lightcones. 

We demonstrate below numerically that the curvature 
at the Cauchy horizon x — Xf is finite in Choptuik's 
scalar field critical solution. Furthermore all geodesies 
cross it in finite affine parameter. The spacetime can 
therefore be extended beyond, but this continuation is 
not unique. Mathematically speaking we shall see that 
the solution is not analytic at the Cauchy horizon in the 
limit coming from the past, and so there is no preferred 
analytic continuation to the future. (The curvature is 
only C° from the past.) The family of continuations 
in which the curvature is C° across the Cauchy horizon 
and which are DSS is parameterized by one free periodic 
function U e (r) . Physically speaking this function can be 
interpreted as data on the naked singularity which de- 
termine the continuation, in addition to the null data on 
the Cauchy horizon. 

We now discuss the possible continuations that are al- 
lowed kinematically. 

In the simplest case, the £-lines to the future of the 
Cauchy horizon are all timelike (A < 0) until a timelike 
singularity F — is reached. In our classification this is a 
dynamical singularity, while the kinematical singularity 
is a single point. This conformal structure is shown in 
Fig. [3 If m ~ r 3 , or equivalently fi ~ F 2 the 
conformal structure is the same, but with the singularity 
replaced by a regular center. 
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FIG. 1: (Uncompactified) conformal diagram of the critical 
solution up to the Cauchy horizon showing DSS-adapted co- 
ordinates. DSS lines are shown continuous. Lines of constant 
r are shown dashed (we have assumed e A = 2 while in reality 
e A ~ 31). Note that the numerical domain is < r < A 
(shaded), with periodic boundary conditions. We have illus- 
trated the three coordinate patches we use for numerical work: 
in the past patch between x c and x v r-lines are spacelike. In 
the outer patch between x p and x; r-lines are timelike. In 
the future patch beyond Xf they are null. The three patches 
together cover the entire spacetime without overlapping, and 
the coordinates x and r are continuous at the interfaces. 



central 
singularity 



regular 
center 




FIG. 2: An extension with a timelike central singularity. In- 
stead of the central singularity the solution could also have a 
regular center. This diagram has two fans at x p and xj. 



FIG. 3: A hypothetical extension with a null singularity. This 
diagram has two fans at x v and xf and a splash at Xh- Note 
that this extension is not realized as a DSS continuation of 
the Choptuik solution. 



Nolan [f| has drawn spacetime diagrams for spherically 
symmetric CSS spacetimes in which the point singular- 
ity at the origin of the Cauchy horizon is only the start- 
ing point of an ingoing central null singularity. In our 
classification this is an extended kincmatical singularity, 
which requires the existence of at least one splash, that 
is a line where A = again. Fig. [3] shows the simplest 
generic possibility, in which the splash is followed by a 
spacelike dynamical singularity, which covers part of the 
naked null singularity. This spacetime structure can ac- 
tually be realized in spherical CSS scalar field solutions 
0- The same, and more exotic structures, can also be 
realized in spherical CSS perfect fluid solutions |g. 

Contrary to our expectations, we have found that none 
of these exotic possibilities are realized as spherical DSS 
continuations of the Choptuik solution. There is a unique 
DSS continuation with a regular center, and all other DSS 
continuations have a timelike singularity. In hindsight, 
the reason appears to be that the null data on the Cauchy 
horizon are extremely weak. With stronger data (not 
associated with the Choptuik spacetime) we find different 
kinds of continuations. 



III. COORDINATES AND FIELD EQUATIONS 

In spherical symmetry, there are four algebraically in- 
dependent Einstein equations, which can be taken to be 
G TT , G rx , G xx and Gee- The fourth of these is a combi- 
nation of derivatives of the first three and can therefore 
be disregarded. The first three equations contain first 
and second derivatives of F, but only first derivatives of 
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A, B, and C. In the following we assume that F is a given 
function of x and r. The Einstein equations then deter- 
mine three independent linear combinations of the six 
first derivatives of A, B and C . A non-degenerate basis 
of this three-dimensional space is n tX , fi^ T and V^W. 
(With F(t, x) given, this last term contains only first 
derivatives of A, B and C.) 

Here we investigate massless real scalar field matter. 
The scalar field <fi obeys the wave equation 

V a V a ^ = 0, (6) 

and the Einstein equations can be written as 



R ab = %-KGV a cf)\J b 



(7) 



The wave equation, when written as a second order PDE, 
will in general contain both first derivatives of all four 
metric coefficients. However, by writing the wave equa- 
tion in a geometrically defined first-order form, all met- 
ric derivatives except one can be eliminated. Define two 
null derivative operators V u and V u with the usual con- 
vention that both point towards the future and that V u 
points inwards while V„ points outwards. Normalize 
them by imposing that V„r = — r and V ' v r = r. Wc 
define the null derivatives of the scalar field as U = 
\j2nG V u (j) and V = \j2itG V.„0. The massless wave 
equation in spherical symmetry can then be written as 



v u v- 



V 

u 



PU 
PV 



with the scalar P = rV 2 r/(Vr) 2 



(8) 
(9) 

1. Using the 

Einstein equations, the curvature can be given in terms 
of U and V: 

Run = 4U 2 , R vv = W 2 , R uv = AUV, (10) 
(other components of the Ricci tensor vanish) and 

-4 



R 



-UV. 



(11) 



The most general form of the scalar field compatible 
with DSS is 



(t, x) — if}{r, x) + KT 



(12) 



where ip(r + A, x) — tp(T,x) and k is a global constant. 
In the Choptuik solution, k = empirically, so that <f> is 
periodic in r. This means that U and V have zero average 
in t. Moreover, U and V in the Choptuik solution obey 
U(t + A/2, x) = -U(t, x) and so for V. (This of course 
implies zero average.) As a consequence /z(r + A/2, x) = 
fi(r, x) and so for other suitable metric fields. We shall 
assume these extra symmetries in our numerical work, 
but all our analytical expressions remain valid if these 
symmetries are dropped. 

We now describe three coordinate patches that cover 
the critical solution. We demand that both the past and 
the future lightcones of the singularity occur at lines of 



constant x. This makes it easier to impose regularity 
at the center and past lightcone, and to investigate the 
behavior at the future lightcone (Cauchy horizon). It 
also allows us to match the coordinate patches without 
overlap. Subject to these requirements, we have tried to 
make the field equations in each patch as simple as pos- 
sible. Based on the three patches, it is straightforward to 
construct a single smooth coordinate system covering the 
whole spacetime (see Section IVH . but using it from the 
beginning would unnecessarily complicate our numerical 
work. We summarize a number of coordinate systems 
for spherically symmetric CSS and DSS spacctimes, and 
their advantages and disadvantages, in the Appendix. 



A. Past patch 

On the past patch, which extends from the regular 
center to the past lightcone, we write the metric in terms 
of two free functions /(r, x) and a(r, x) as 



A = a 2 (x 2 -.n 



B = -XCL 



C = a 2 



F = x, 
(13) 

for x > x c = 0. Here a = (1 — 2m/ r) x l 2 is the scalar 
we defined above. These coordinates can be derived from 
the Schwarzschild-type metric 



ds 2 = -a 2 dt 2 



1 j 2 

i dr 



r 2 dn 2 



(14) 



through the coordinate change x — r/(— t) and r = 
— ln(— t), and defining / = a/a. The remaining gauge 
freedom of time relabelling is fixed by imposing that the 
past lightcone is on a constant x = x p line. Note that 
we do not follow the convention of that a = 1 at the 
center, which in |5j forced us to introduce an additional 
free function of time in the definition of x. The Einstein 
equations are 



f,x 
(«~ 2 ),x 



(a 2 ~ 1)/ 

x 

l-(l + U 2 + V 2 )g- 2 

X 

{f + x)U 2 - (f-x)V 2 



1 



and the matter equations are 

/[(l - a 2 )U + V] -xU, T 



U x = 
V x = 



x(f + x) 
/[(l -a 2 )V + U] +xV T 
x(f - x) 



(15) 
(16) 

cr 2 - 1(17) 



(18) 
(19) 



At the regular center x — we impose elementary 
flatness, that is the absence of a conical singularity. In 
order to do this, we define H — (V + U)/(2x) and ^ = 
(V — U)/(2x 2 ) and impose that both are regular even 
functions of x at x — 0. At the past lightcone we have 
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/— x = 0, which by our gauge choice means that / = x p 
there. We also impose the physical regularity condition 

V T - (1 - a 2 )V~ U = (20) 

on the past lightcone. The conditions of DSS, regular- 
ity at the center and regularity at the past lightcone se- 
lect a solution. The equations on the past patch are 
form-invariant under the linear coordinate transforma- 
tion x — > cx, f — ► cf, and a, U and V unchanged. In the 
numerical results presented here we have set x p = 1 on 
the past patch. Note that the regularity condition 12U|) 
is coordinate- independent, as U, V and a and r are all 
scalars and r = — lnr on the lightcone. 



B. Outer patch 

On the outer patch, which extends from the past to the 
future lightcone, we write the metric in terms of a(r, x), 
b(r,x), and the auxiliary function £(r) as 

A = a 2 {l-b 2 ), B = ab£„ C = -£ 2 , F = 1, 

( 21 ) 

where £ > is a function of r only, a > is the scalar 
defined above. We fix the remaining gauge freedom by 
imposing that the past lightcone b = — 1 occurs at x = x p , 
and the future lightcone b = 1 at x = xt. Putting both 
lightcones on a line of constant x requires £ (r) to be non- 
constant. There is no outer coordinate patch that does 
not have at least one function like £(t). As t = — lnr 
everywhere on the outer patch, r is continuous between 
the past and outer patches. 
The metric equations are 



-3 + a 2 + U 2 {l-b) + V 2 {l + b) £' 



2a 



U 2 -V 2 



£ 



(23) 
(24) 

(25) 
(26) 



{a- 2 ) r = [l + U 2 (l-b) + V 2 (l + b)] a- 2 - 
and the matter equations are 

Up _ (1 -a 2 )U + V + U iT 

~t~ ^0~b) s 

V x (l-a 2 )V + U + V T 
T ~ a(l + b) ■ 

Note that the a tT constraint equation is again linear in 
a~ 2 . The equations on the outer patch are form-invariant 
under the linear coordinate transformation x — > cx + d, 
£ — > £/c, and a, 6, ?7 and V unchanged. In the numerical 
results presented here we have set x p — — 1 and Xf = 1 
on the outer patch. 



C. Singular behavior at the Cauchy horizon 

We shall see that at the Cauchy horizon the solution 
is mildly singular. Naive finite differencing breaks down 



there. Instead we expand the generic solution around the 
Cauchy horizon in terms of two free periodic functions 
Vo and U e , and match this expansion to the numerical 
evolution at a small finite distance to the past of the 
Cauchy horizon. Before describing the full expansion, we 
focus on the origin of the singular behavior. 

Equation (|25(l becomes singular at the future lightcone 
because the denominator of the right hand side vanishes 
there. In contrast to the past lightcone, we do not have 
any freedom left to enforce the vanishing of the numer- 
ator as well. Therefore we expect the solution to have 
some kind of singularity at Xf . We rewrite (|25|l 



DU X - U T - (1 - a 2 )U = V, 

where we have defined the metric function 

a(l-b) 



D 



£ 



(27) 



(28) 



D is positive on the outer patch and vanishes on the 
future light cone. The characteristics x(t) of Eq. l(2T|) 
are given by 



dx(r) 
dr 



-D (t,x(t)) 



(29) 



Recall that we impose the gauge condition that the CH 
is at x = Xf, or b(r,Xf) — 1. We define the shorthand 
y = x- x f . 

We now make one fundamental assumption, namely 
that the spacetime admits regular null data on the 
Cauchy horizon y — 0. This assumption is clearly nec- 
essary if we want to continue the spacetime through the 
CH, but here we make it simply because we have not 
been able to find a more general ansatz. We shall see 
later that it is sufficiently general to be matched to the 
critical solution that we have obtained numerically on the 
past patch. 

We therefore assume that V is continuous, or 



V(t,x) = V (r)+o(y ). 



(30) 



By substituting this into Eq. 124|) in the limit b = 1 we 
find that a is also continuous and oo(t) obeys 



(ao 2 )'-(l + 2^ 2 )a - 2 + l = 0. 



(31) 



Because we impose periodic boundary conditions in r this 
ODE has a unique solution. The physical significance of 
this is that the null data Vq determine the geometry of the 
CH. Similarly, from (|27|l in the limit D = we find that 
U must be continuous, and Uq(t) is the unique periodic 
solution of 



U' + {1- a 2 )U +V = 0. 



(32) 



This condition follows from the assumption of DSS. Fi- 
nally, from (J2EI, 123 and (J2H| we find that D is once 
differentiable, so that D(t 7 x) — y D\{t) + o(y), and D\ 
is given by 



£ 



2V n J 



(33) 
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At this point we introduce more shorthand notation. If 
/(r) is any periodic function (with period A), let / be its 
average value, and let /(r) = /(r) — / be its oscillatory 
part. Let f f be the definite integral f(r') dr' where 

To is chosen so that J f has vanishing average. 

We can now integrate l|29(l for the U -characteristics to 
leading order and obtain 



log \y{r)\ + D\T +/£>! + o(y°) = const. 



(34) 



We see that on a characteristic r ^ — oo as |y| ^ 0. 
Because U (r, x) is periodic in r with period A, an infinite 
number of oscillations in y at constant r pile up at the 
Cauchy horizon y — 0. We can solve <|27[) to leading order 
by the method of characteristics. The general solution is 

U(t, y) = U (t) + \y\ e Ue(T) U e (t) + o(|y| e ) (35) 

where U e (f ) is an arbitrary periodic function with period 
A and 



T 

H(t) 
K 



1-al 



exp 



t + H(t)+ Khx \y\, 

Di J 
1 



Rewriting l|31|l as 

a2-(l + 2^ 2 ) = 2(lna )', 



(36) 

(37) 
(38) 
(39) 

(40) 



(41) 



we see that a 2 , = 1 + 2Vg 2 , and so we can express e and 
K as 



1 - 2V 2 



K 



-(1 



(42) 



Our initial assumption that U and V are continuous at 
the light cone is therefore justified either if < Vq < 1/2, 
so that e > 0, or if U e (r) = 0. We shall show numerically 
that e is small but positive on the CH. In this case U 
and V are just C°, and the scalar field is therefore C 1 . 
In spherical symmetry the Riemann tensor is determined 
completely by the Ricci tensor, which in turn is quadratic 
in the partial derivatives of 0, see (J7J). The curvature is 
therefore quadratic in U and V and so is C°. 

A similar analysis, with U and V interchanged, applies 
to the past lightcone. In the notation we have introduced 
here we can describe the past lightcone by saying that 
e < there (because the null data Uq on the past light- 
cone are large, Uq > 1/2) but that the free coefficient V e 
vanishes identically (because we have imposed analyticity 
as a boundary condition). 



D. Expansion near the Cauchy horizon 

We cannot apply Fuchsian techniques to our system of 
equations because they require the simultaneous vanish- 
ing on y = of the coefficients of U x and U T in equation 

However, the form (|35|l of the leading terms in U sug- 
gests that the full non-linear solution can be written as 
a regular part, containing only integer powers of y, plus 
a singular part which contains powers of |y| e . We can in 
fact construct a formal solution near the Cauchy horizon 
in the form of an asymptotic double series: 

OO OO femax(n) 

/(^) = E» n /«(r) + E E \v\ n+k£ fn + Ur,x), 

(43) 



n— 1 n— 1 k— 

where / stands for U, V, a and b, and 



f n+ke (T,x) 



,ax(n,fe) 

V f {€> 



(r)B ke (f)- (44) 



Here e and f are the quantities defined above in terms 
of Vq. The expansion depends on the two free periodic 
functions Vq(t) and U e (f). By function counting we can 
therefore match any DSS solution to this expansion. 

The first non-integer term appears in each variable at 
the following orders: 

U(t,x) = U (T) + \y\m e (T)U e (r) 

+yU 1 (r)+0{\y\ 1+e ), (45) 

V(t,x) = Vfc(r)+tfVi(T) + 0(|tf| 1+e ), (46) 
a(T,x) = ao(r)+yai{r) + 0(\y\ 1+e ), (47) 
b(r,x) = l + yb 1 (T) + y 2 b 2 (T) + 0(\y\ 2+ '). (48) 

In the previous section we obtained the coefficients of 
expansions (|45H48|) up to 0(\y\ e ). Stopping there, the 
first order we neglect is 0(y). This truncation already 
depends on both free functions Vq and U t and shows the 
singular behavior. It is also a sensible truncation numeri- 
cally because e turns out to be very small in the Choptuik 
solution. Going further, for the same reason there would 
be no point in including 0(y) terms without also includ- 
ing all 0(|y| 1+fce ) terms. It turns out that we need to 
go to 0(|y| 1+3c ). We have used the expansion to that 
order to check convergence. The expressions are given in 
Appendix 



E. Future patch 

Our analysis of the possible continuations of the criti- 
cal solution in Section[H]has shown that we can cover the 
entire future of the Cauchy horizon in a single patch if we 
make r an ingoing null coordinate. This means setting 
C = and B < 0. In order to put the center r = at 
a known coordinate location, we also set F = —x. We 
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choose x < here so that x increases as we extend the 
spacetime away from the Cauchy horizon. We parame- 
terize this metric in terms of the scalar a and a coefficient 
/ (not the same as / in the past patch): 

A = -4a 2 f(f + x), B = 2a 2 f, C = 0, F = —x. 

(49) 

Regularity of the metric requires a > and / > 0. The 
field equations are 



(a~ 2 X* 
U, x 



{o 2 - 1)/ 

X 


(50) 


1 - (f + 2U 2 )a- 2 


(51) 


) 

X 


r 2fv 2 -2(f + x)u 2 | x 

X 


a~ 2 - 1,(52) 


f[{\- a 2 )U + V] -xU T 


(53) 


x(f + x) 


U+(l- a 2 )V 


(54) 



At the Cauchy horizon x = Xf < we impose the coor- 
dinate condition / + x = 0. The equations are invariant 
under x — ► cx, / — * cf. We use this to set Xf = — 1 
on this patch. As a consequence r = — In r on the light- 
cone, and so t is continuous between the outer and future 
patches. 

We find an asymptotic expansion around the Cauchy 
horizon in terms of two free functions Vq and U e . Vq is 
given by the null data on the Cauchy horizon, and so is 
the same as on the outer patch, e is therefore the same 
on both sides. Uq obeys the same ODE, (|32[1 . as on the 
outer patch, and so is the same function. There is no 
need, however, to make U e the same on both sides, as we 
would not gain any differentiability by doing so. Instead 
we consider U e as free "data on the naked singularity", 
and we shall find experimentally how the global structure 
of the spacetime is influenced by this choice. 



IV. NUMERICAL CONSTRUCTION OF THE 
CHOPTUIK SPACETIME UP TO THE CAUCHY 
HORIZON 

We have carried out a brand new numerical computa- 
tion of the Choptuik spacetime, improving the precision 
of our previous calculations by roughly four orders 
of magnitude. This was mainly needed to assert with- 
out doubt that e is different from zero, even though it is 
extremely small. 

Essentially, our new scheme uses shooting methods on 
the x axis, instead of relaxation methods. This allows us 
to improve the treatment of the regular singular points 
of the equations (the center and the lightcones) by using 
Taylor expansions. We still work with pseudo-spectral 
Fourier techniques in r because the solution is periodic. 
We shall see however that the particular structure of the 
Fourier transform of the Choptuik spacetime poses an 



unexpected problem when combining Taylor and Fourier 
expansions. 

In this Section we explain the numerical scheme and 
present the results for the first two patches, with special 
emphasis on convergence properties and error analysis. 



A. Numerics 

1. Pseudo-spectral decomposition 

For definiteness, let us suppose we work on the past 
patch. Following |5j we discretize our A-periodic fields 
Z(t,x), where Z stands for any of the set {a, f, U, V}, 
using N equidistant points in one period: 



N—i 



Z n (x) = Z(-A,x) = Zk( 



(55) 



for n — 0, N — 1. In this way we transform our 1+1 
PDE problem for Z(t, x) into an ODE problem for the 
modes Zk{x). The essential idea of pseudo-spectral meth- 
ods is to carry out algebraic operations pointwise on the 
Z n and T-differentiation/integration on the Zk, switch- 
ing from one to the other with a Fast Fourier Transform 
algorithm. 

The main drawback of the method is the aliasing prob- 
lem: pointwise products of fields (nonlinearities) gener- 
ate high frequency modes which cannot be sampled with 
only iV points. We partially solve that problem by dou- 
bling the Zk (padding with zeros) before going to the 
Z n , carrying out the necessary algebraic operations on 
the doubled Z n and going back to the Zk, then halving 
the Zk and thus throwing away high frequency noise. We 
have tried other possibilities, such as padding with 3N 
or 7N instead of N zero Fourier components, or extrap- 
olating the Fourier coefficients using the observed fact 
that high frequency modes have a simple exponential de- 
pendence on frequency (see below), but the results are 
not improved. Aliasing can only be reduced by going to 
higher N. From a numerical point of view, we are only 
safe from aliasing when the amplitude of the modes we 
are cutting off is below machine precision. 

Because all our fields are real Z^ = Z^-h- Further- 
more, the metric fields a and / are even [in the sense 
a(r + A/2, x) — a(r, x)] and therefore their Fourier trans- 
form only contains even fc modes, while the matter fields 
U and V are odd [in the sense U(r + A/2, x) = —U(r,x)] 
and therefore their Fourier transform only contains odd 
k modes. Taking this symmetry into account an even 
or odd field Z sampled with N points per period is en- 
coded by N/4 independent non-zero complex modes. As 
we have 4 independent variables a,f,U, V, the ODE sys- 
tem we solve comprises N complex or 2N real variables. 
In our calculations we have used N = 32, 64, 128, 256 and 
512. Previous investigations used N = 64. See Appendix 
C of Ref. p| for a complete discussion of our Fourier 
pseudo-spectral method. 
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2. Shooting to fitting points 

We cannot cross the lightcones during the integration 
of the ODE in the x-axis because they are regular singular 
points of the equations. Therefore we perform consecu- 
tive shooting calculations on the past, outer and future 
patches, in this order because we need information from 
the first to shoot the second and from the second to shoot 
the third. The issue of error propagation becomes very 
important. 

Again, we describe the past patch for definitencss. 
Given the equations H15I - I18J) and free data /(r, x c ) = 
/ c (t), ^(t, x c ) = ^ c (t) at the center, we calculate the so- 
lution at a;ieft slightly larger than x c using a second-order 
power expansion [leaving errors of order (:ci ft — a; c ) 4 ]. 
From these data we integrate the ODE system forward 
in x using finite differencing, up to x m id- In the same 
way, given free data U(t, x p ) = U p (t) and the gauge con- 
dition f(r,x p ) — x p , we calculate the solution at x r ight 
slightly smaller than x p [this time with errors of order 
(x p — Xright) 3 ] and integrate backward in x up to the 
same £ m id- Finally we use Newton's method to look for 
the free data which brings the mismatch between both in- 
tegrated solutions at x m ;d down to a minimum, typically 
of order 10~ 13 . (This is the machine precision of 10~ 16 
reduced by a factor 10 due to the calculations at each 
step and a factor 100 from the ODE-integration along 
sa 10 4 points.) 

We use a grid that becomes logarithmic near the regu- 
lar singular points, with maximum stepsize /i m ax at some 
intermediate position. That grid is constructed using the 
transformation 



l + e z 



(56) 



from a grid of equidistant points in z between the val- 
ues z\ e { t and bright corresponding to xi oft and 

bright re 

spectively. Near the two endpoints we have x — x c ~ 
(x p — x c )e z and x p — x ~ (x p — x c )e~ z . We integrate 
on a fixed grid in x rather than using a variable stepsize 
method in order to check for convergence with h max . This 
gives us a good estimate of the underlying discretization 
error. 

Concerning the ODE-integrator, we have tried several 
Runge-Kutta methods, both explicit and implicit (Gauss- 
Legendre), with different convergence orders In gen- 
eral, implicit methods are better suited to our problem 
than explicit methods, particularly for high N, because 
high frequencies make the problem stiffer. We choose 
an IRK2 method (implicit, 2 stages, order 4), which is a 
compromise between the accuracy of a high order method 
and the clarity of convergence of a low order method. 
Our implicit schemes arc implemented by iteration un- 
til ^-differences between successive iterations converge 
below 10 -15 . We cannot get closer to actual machine 
precision (10 -16 ) due to the accumulation of roundoff er- 
ror. An implicit step typically takes 5 to 15 iterations to 
converge. 
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FIG. 4: Best free data on the singular points: / c (t) (con- 
tinuous line), \& c (r)/300 (dotted line) and U p (t) (dash-dot 
line) . 



B. Past patch 

The natural choice for the coordinate of the center is 
x c = 0, but there is no preferred value for the past light- 



cone coordinate and we choose x 
parameters 



p 



1. With the set of 



N 

bright 



256, 
0.001, 
0.01, 
0.9999, 
7- 10~ 4 , 



(57) 



our Newton's method converges to the free data plotted 
in Fig. 01 (see also table P), with a value 



A = 3.445452402(3), 



(58) 



which improves the precision of our previous result 
3.4453(5) more than four orders of magnitude. The met- 
ric and matter fields integrated from those free data are 
given in Fig. 

The error bars in Q58J1 and tableware estimated from 
the convergence properties of the code, which we now ex- 
plain in detail. In general, we have observed that we can 
reach higher relative precision in the metric fields than 
in the matter fields, because the former are essentially 
integrals of the latter. Figs. El and |S1 contain all the con- 
vergence information for the past patch, but in order to 
properly discuss convergence issues, we first need to talk 
about an important feature of the Choptuik spacetime. 

Fig-Hshows the field ^'(t, x) together with a logio plot 
of the modulus of its Fourier transform in r. There is a 
clear difference in behaviour between the regions above 
and below x ~ 0.2. (This difference is present in all our 
fields, but it is particularly important in ty, as we will 
see.) Near the center the function has large r-derivatives 
which require many modes in the Fourier expansion to 
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TABLE I: The first 8 nontrivial Fourier modes of the free data. (The error in the last digits is shown in brackets.) Note that 
some of them have a relative precision better than 10~ 8 . 



Re f 2k (xc 



Im f 2k jxc) 



Re <j> 2 fc+i(zc) Imj>gfc+i(gc) ReU 2k+1 (x p ) 



Im U2k+i(x p ) 



0.2071909728(5) 
by def. 
-0.02370998706(19) 
0.01347536638(18) 
-0.001117368391(5) 
-0.00432309030(6) 
0.00345031858(11) 
-5.3677559(20) 10" 4 





0.0649057078(3) 

-0.01438139603(19) 

-0.00645699456(5) 

0.00883071824(10) 

-0.00355712461(5) 

-0.00117027642(11) 

0.00236648443(8) 



0.788624247(31) 

11.52753960(12) 

-9.12055613(23) 

-1.5735907(5) 

8.2142645(16) 

-5.8935621(21) 

-0.6075901(6) 

4.3606105(27) 



-11.6194821(5) 

4.2699131(6) 

7.5700655(6) 

-10.3634715(8) 

3.3866919(5) 

4.3873917(11) 

-5.9471050(31) 

2.0264097(27) 



0.2962634507(4) 
-0.00901909339(14) 
-0.002037853110(17) 
2.55305777(10) 10~ 4 
1.12390710(11) 10~ 5 
-4.9301780(6) 10~ 6 
1.884812(9) 10 -7 
7.04459(5) 10~ 8 



0.0905094329(7) 
-0.02200156878(7) 
0.00159029448(4) 
1.73238389(4) 10~ 4 
-3.6797294(3) 10~ 5 
8.2486(15) 10" 9 
6.14734(5) 10" 7 
-4.73630(8) 10~ 8 





FIG. 5: Past patch fields on a single A-period. 



be resolved (see also Fig. . Far from the center those 
derivatives are much smaller and just a few modes are 
enough to achieve high resolution results. This is re- 
flected in a very slow decay of the Fourier transform near 
x ~ x c and a fast decay for x above 0.2, although both 
are exponential decays. It is important to emphasize 
that this exponential decay is the best numerical evi- 
dence we have in support of the analytical character of 
the Choptuik spacetime, given the absence of a math- 
ematical proof of existence of an analytical solution to 



which our numerical spacetime should be converging. 

We have observed as well that the phases of the high 
frequency modes Z k tend to a linear dependence on k at 
a given point x. Therefore the high frequency behaviour 
of our fields Z is similar to that of the function 



E 



- a I k\ 



sinha 



cosh a — cos ^g 2 - ' 



(59) 



which has periodic sharp peaks of width aA/(47r) and 
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FIG. 6: Convergence figures for the past patch. The three columns correspond to f c , ^ c and U p , respectively. The six rows are 
organized in three pairs, describing convergence with respect to /imax, xieft and bright respectively, when these three parameters 
are halved two or three times. The first row of each pair shows consecutive differences of fields, rescaled by factors 4, 16 and 
5.2, respectively, so that they coincide when converging with orders 2, 4 and 2.4, respectively. (Convergence with respect to 
bright is slower than the expected order 3, in particular that of the very low frequency modes of f c .) The second row of each 
pair shows a log 10 of the power spectrum of those consecutive differences (without rescaling), to show the different behaviour 
of the Fourier modes. Convergence of the high frequency modes of 9 C is worse than that of their low frequency counterparts, 
as explained in the text. 
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FIG. 7: \1/(t, x) and its power spectrum | ^ 2 k+i(x) | . Note the difference in behaviour between the regions below and above 
x~0.2. 



height 2/a for small a. In the Choptuik spacetime the 
decay exponent a ranges from 0.28 for f c at the center to 
2.04 for Up at the past lightcone. We have not found an 
explanation for this behavior, but it seems to be of dy- 
namical origin: Arbitrary high-frequency perturbations 
of the correct free data at the center decay towards larger 
x, and high-frequency perturbations at the past lightcone 
grow when integrated towards the center, probably due 
to 1 jr factors in the equations of motion. 

From a numerical point of view, this means that alias- 
ing problems will be important near the center x = 0. 
This forces us to use a large TV (TV — 512 in Fig. 01. But 
then most of the modes become essentially random for 
x > 0.2 because they are well below the error thresh- 
old given by roundoff error, estimated in 10~ 14 in rela- 
tive terms. This threshold generates the flat plateau in 
the right panel of Fig. The main sources of round- 
off error are threefold: ODE integration along x, Fourier 
transforms, and inversion of a very stiff matrix in New- 
ton's method. Particularly important are the errors in 
the modes of U p above k — 15 because they propagate 
inwards and get amplified, giving errors of relative or- 
der 10 -9 in the matter fields, mainly in $> c . (See Fig.[7| 
again.) This fact sets the limit of the maximum accuracy 
that we can get in the results, with just double precision 
numerics and using our code. We could use quadruple 
precision to improve the precision in U p but the calcu- 
lations would become too slow. Alternatively, we could 
force the vanishing of those modes that we believe must 
vanish, but we do not want to assume anything at all 
about the result in advance. 

We now check convergence with respect to the numer- 
ical parameters. As one would expect, the final solution 
is completely insensitive to the choice of intermediate fit- 
ting point a; m idj although the convergence of Newton's 
method is faster when using smaller values because the 
mismatches are larger. We choose x m id = 0.01. 

Several tests in simpler problems show that the ODE 
integrator in x is perfectly fourth-order convergent. The 



first two rows of figure also show this fact, even though 
roundoff errors slightly blur the point. Note that the 
modes of U p (t), after the first 10, do not converge be- 
cause they are already below our error threshold (on the 
plateau in Fig. [7J| . Note also that high- frequency modes 
of ^ c do not converge for such small values of /i max - They 
do converge for larger values. 

Convergence with TV is exponential as expected. Fig.|5| 
shows the power spectra of the free data for TV — 
32, 64, 128, 256, 512. It clearly shows how much the re- 
sults are improved by doubling TV and how the plateau 
goes down each time, until TV — 256, when errors in U p 
hit our error threshold. The data for TV = 512 shows 
that we cannot improve the results any further because 
we cannot decrease the errors in U p . f c is then perfectly 
resolved down to machine precision, but high-frequencies 
of the matter variables cannot be improved near the cen- 
ter. It is clear that using TV = 256 or TV = 512 we are 
not affected by aliasing errors. 

Finally, convergence with xi c ft and alight depends on 
the value of TV: As we said, we calculate the fields at x\ e f t 
and a; r ight using Taylor expansions: 

/(t, x) = f c (r) + f 2 (r)(x - x c f +0[{x- x c ) A ] (,60) 
U{t,x) = U v (t) + Ui{t)(x-x p ) 

+U 2 (t)(x - xp) 2 + O [(x - x p f] , (61) 

and so for the other fields. Therefore we expect fourth 
order convergence with respect to xi e ft and third order 
with respect to alight- The coefficients of the Taylor ex- 
pansions are obtained as nonlinear combinations of the 
free data and their r-derivatives. The latter are calcu- 
lated multiplying the Fourier modes by ik, which am- 
plifies the high frequency modes. Fig.[§]shows the Fourier 
transforms of $ c ,$2j^4 obtained for different values of 
TV. For low TV our estimations of these coefficients are 
bad and we do not see the expected convergence order 
in the expansions with x\ e ft and £ r i g ht, both because we 
are cutting off too soon in frequency, leaving out modes 
which are important (see for instance the case of ^2 with 
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10 +2 




mode 



fore we can estimate the error of the base run (|57(l with 
respect to any of the parameters of the code as the differ- 
ence between the base run and another run with a refined 
value of that parameter. Those are precisely the contin- 
uous lines in Fig. [5] Truncation errors from space and 
time discretization could be reduced to machine preci- 
sion. However, errors from the expansions at the singu- 
lar points cannot be brought down with our code below 
a limit which we estimate (assuming no systematic er- 
ror) between 10 -8 and 10 -9 in relative terms. Therefore 
there is no point in reducing the other sources of error 
below that limit, and the choice of parameters for our 
base run l|57l) reflects this. 



FIG. 8: Power spectra of the free data for 8, 16, 32, 64 and 
128 modes. Note the very different behaviors of the functions 
f c and ^ c at the center from the function U p at the past 
lightcone. 
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FIG. 9: Power spectra of the functions $ c , ^2 and *l/4 for 
8, 16, 32, 64 and 128 modes. Aliasing problems show up as 
unphysically growing tails for the highest frequencies in each 
case. Working with only 32 modes, this introduces errors in 
^2 which are of the same order of the amplitude of the most 
important modes. We cannot avoid going to at least 64 modes 
(that is, N = 256). ^4 is presented for illustration purposes 
only; we do not use it in the code. 



N = 64), and because of aliasing errors, which gives us 
wrong estimations of the modes that we are including 
(see the unphysical tails at the end of the functions). 
The same phenomenon happens on the past lightcone, 
but its effect is not so important. That is the main rea- 
son why we need at least N — 256 to get good results. 
However, for higher N we do see clear convergence with 
the expected orders (with the exception of the very low k 
modes, whose convergence with respect to alight is slower 
due to accumulation of high-frequency errors in U p ). This 
is shown in Fig. 

We conclude that the code converges in the expected 
manner with respect to all numerical parameters. There- 



C. Outer patch 

In the entire outer patch the fields are much smoother 
in r (as they are already on the past patch for x > 0.2). 
Therefore we only need a few Fourier modes to get good 
precision. In fact N = 64 is enough to reach the maxi- 
mum accuracy given by propagation of the errors in the 
past patch. There is no reason to increase N further. 

With x p = — 1 and Xf = 1 we choose these parameters 
for the numerical evolution: 

N = 64, 
xieft = -0.9999, 

X m id = -0.9, (62) 

bright = 0.9999, 
h max = 0.001, 

Now the free data are the metric function £(r) and the 
matter functions V/(r) = V(r,Xf) and U e (r) at the fu- 
ture lightcone. (Note that the function Vf was called Vo 
in Subsection IIII (JI )The results are given in table UTI and 
Fig. ^3 with a final value 

e= 1.4710439(8)- 10" 6 , (63) 

clearly different from zero. We firmly believe that this 
is neither a numerical artifact nor a consequence of our 
expansion around the future lightcone. After showing 
convergence of the code in the outer patch, we dedicate 
most of this subsection to supporting this claim. 

The integrated functions are shown in Fig. ^2 The 
runaway of characteristics is not apparent in this figure, 
because all but the first oscillations are piled up in the 
region between x — 0.95 and x = 1. Fig. 1121 shows U 
using a logarithmic i-axis. 

We first analyze the issue of error propagation from 
the past patch to the outer patch. We need to find out 
how small variations in A and U p change the results of 
the outer patch shooting. Assuming that the former are 
small, we calculate first derivatives of the latter. The 
variations of the outer patch free data with respect to A 
are: 

il = 5. 4 .l0- 6 , l^-Hog = 0.0020, (64) 
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TABLE II: The first 8 nontrivial Fourier modes of the outer patch free data. (The error in the last digit is shown in brackets.) 
Note the very different relative precisions achieved in £ and U t . The former is insensitive to changes in the parameters of the 
code, but not the latter. 



Re £2 



Im £ 2 



Re V 2 k+i(x f ) 



lmV 2k+1 (x f ) ReU e2k+1 ImU e2k+1 



1.322045988(6) 
-0.00492853319(12) 
-1.1237092(5) TO" 5 
1.99072448(4)T0" 5 
-1.078319002(12)T0"' 
-1.032893598(26) -10" 
1.74888526(8) TO" 8 
2.39356(8) TO" 11 





-0.00430580184(14) 
2.89549492(5)T0~ 4 
-6.20831744(26) TO" 
-1.46908361(3) TO" 6 
1.436706184(10) -10^ 
5.7084263(10) TO -9 
-2.0053466(8) TO -9 



-5.177664(12)T0 
-1.5(3)T0- 10 
-1.(4) TO" 12 
-1.(5) TO" 13 
1.(30) TO" 14 
-0.(21) TO" 14 
1.(25) TO" 14 
0.(3) TO" 13 



■ 4 3.157186(5)T0" 
-2.0(7) TO" 10 
2.(15) TO -13 
0.(31) TO" 14 
-1.(6) TO" 13 
0.(4) TO" 13 
0.(4) TO" 13 
0.(3) TO" 14 



-0.03844(5) 
0.003883(8) 
2.895(8)T0" 4 
1.748(5) TO" 5 
8.957(26) TO" 7 
4.126(10) TO" 8 
1.826(4) TO -9 
7.79(12) TO" 11 



■0.250325(8) 
-0.0121478(20) 
-5.853(3)T0" 4 
-2.5842(28) TO" 5 
-1.0219(20) TO" 6 
-3.920(10) TO" 8 
-1.454(4) TO" 9 
-5.21(26) TO" 11 
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FIG. 10: Free data in the outer patch: 20V/ (continuous 
line), Ue/20 (dotted line) and £— 1.3 (dash-dotted line). Note 
that they are all quite smooth. 



\m\\c 

8A 



= 0.15, 



6A 



= 0.028. 



(65) 



On the other hand, Fig. 1131 shows the maximum vari- 
ations of the free data with respect to changes of the 
Fourier modes of U p . We see that Vf and £ are only 
sensitive to the very low frequency modes of U p , but U e 
changes with every mode. In any case, every derivative is 
small enough: The largest error bars come from the un- 
certainty in A, and then from those of the first two modes 
in Up. The rest of the modes are practically irrelevant 
for error analysis. This sets the maximum accuracy that 
we can achieve on our final results. Assuming quadratic 
error propagation it is 
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It is still very good due to the tendency of small pertur- 
bations to decay when integrating towards larger x, and, 
as we said, achievable already with N = 64. 



We now analyze convergence in the outer patch. Con- 
vergence with h max in the IRK2 method shows per- 
fect fourth order again. Convergence with iri e ft is ap- 
proximately third order as expected because we expand 
around the past lightcone with a second-order Taylor se- 
ries. See Fig. El Finally, we have performed calculations 
expanding around x f using only the order zero terms and 
including the first order terms. 

Vf and £ converge with £ r i g ht to first order when the 
expansion around the CH is truncated at 0(|y| e ), and 
converge to second order when the expansion is truncated 
at 0(\y\ 1+3e ). This is the expected behavior. However, at 
the same time U e converges to first order in both cases 
(see Fig. 1140 . This indicates that adding the terms of 
order 0(\y\ 1+ke ) (with k = 0, 1, 2, 3) to the expansion still 
leaves some 0(\y\ 1+ke ) error. We are confident that this 
is not a simple algebraic mistake in the expansion. We 
note that the excess error in the periodic function U e is 
entirely an error in its overall phase. We therefore suspect 
intuitively that the runaway phase U ~ U c (t — \n\y\) of 
the solution is to blame, but we have not been able to 
formulate this idea consistently. 

In order to show that e is really different from zero, 
we have to analyze the behavior of the function V with 
respect to xi e ft and x r i g ht- The function Vf has a very 
rapidly decaying Fourier spectrum, as shown in Fig. 1151 
As xidt — > x p , all its Fourier modes converge to zero 
except for the first two, and the amplitude of the second 
mode is more than six orders of magnitud below that of 
the first one (see also Table HT}) . 

Fig. 1161 shows V 2 (x) for several evolutions from the 
same x\ e n = — 0.9999 to 15 different values of alight- The 
agreement is very good. This shows that our expansion 
around the future lightcone (including the singular terms 
given in the appendix) captures the behavior of the so- 
lution. 

We conclude that, with errors of order 10" 9 , 

V f (r) = -1.2128648(22)T0" 3 -cos f — + 0.5475726(13) 

(70)' 

Finally, it is interesting to see that going to 0(|y| 1+3e ) 
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FIG. 11: Outer patch fields on a single A-period. 

is not essential for obtaining an accurate result. Table Hill 
compares the results for e using a zeroth-order expansion 
(that is, we only include the terms of orders y° and \y\ e ) 
with those from the first order expansion (which also in- 
cludes orders y and |2/| 1+fcc ). It is clear that they both 
converge to the same number, even though with very dif- 
ferent rates of convergence. 



V. CONTINUATION ACROSS THE CAUCHY 
HORIZON: THE FUTURE PATCH 

A. The continuation with a regular center 

We cannot continue the solution as flat empty 
Minkowski spacetime after the Cauchy horizon because 
we have a small amount of ingoing radiation there, given 
by Vf(r). The simplest continuation to look for is one 
with a regular timelike center, so that the conformal di- 
agram is the same as for Minkowski spacetime. In this 
case both U and V are small on the entire future patch, 



and we can obtain an approximate solution in perturba- 
tion theory around Minkowski space, using the magni- 
tude e 1 / 2 of Vf as the small parameter. To leading order 
in e we obtain the d'Alembert solution on flat spacetime: 



/ = I + 0(e) 
a = l + 0(e) 

U = e 1 / 2 



(71) 
(72) 



F'(f) 



F(f) + G(r) 



+ 0(e 3 / 2 ), (73) 



V 



cl/2 



-G'(t) + (1 + x) 



F(t)+G(t) 



0(e 3/ f34) 



with t = t — ln(l + x). In order to match the null data 
on the Cauchy horizon x = —1, we need — e 1 / 2 G"(r) = 
Vf(r). Recall that V/(r) is given numerically by 170fl . 
If we want to have a regular solution at the center we 
need F(t) — —G(t). The solution is then completely 
determined, and it is clear that a nonlinear solution exists 
of which this is the leading order, and which can be found 
numerically. 

Because the null data Vf on the Cauchy horizon are 
DSS, it appears highly unlikely that there is another con- 
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.999999 



FIG. 12: U in the outer patch using a logarithmic rr-axis. 
The oscillations are clear, but not their slow decay. 



TABLE III: Convergence of 10 6 e with alight using only the 
lowest order terms in the regular expansion plus only the 
0(y e ) in U (first column) and using the full next order as 
well (second column). Convergence is much faster in the sec- 
ond case, as expected, but it is clear that both converge to 
the same number, within our error bars. Recall that digits 
starting from the fifth decimal are not relevant due to propa- 
gation of errors from the past patch. They are shown in the 
second column to make convergence clear. 



bright 


zeroth order 


first order 


0.8976 


238.46 


0.87085 


0.9488 


66.166 


1.36313 


0.9744 


26.743 


1.49818 


0.9872 


9.4840 


1.48105 


0.9936 


1.4699 


1.47058 


0.9968 


0.5058 


1.470305 


0.9984 


1.4386 


1.470949 


0.9992 


1.7639 


1.4710727 


0.9996 


1.5837 


1.4710534 


0.9998 


1.4355 


1.47104333 


0.9999 


1.4332 


1.47104318 


0.99995 


1.4687 


1.471043825 


0.999975 


1.4799 


1.471043941 


0.9999875 


1.4743 


1.471043921 


0.999999 


1.4713 


1.471043911 



tinuation with a regular timelike center that is not DSS. 
We have obtained the (probably) unique regular contin- 
uation by shooting from expansions around the Cauchy 
horizon and a regular center. The free data for the shoot- 
ing algorithm (given by C/ £ at the CH and /, U. x at the 
center) were obtained from the flat spacctime approxi- 
mation. In this case we use an IRK1 integrator and 

N = 16, 
XMt = -0.999, 



bright 



-0.2, 
-0.001, 
2.5 • 10 _ 



(75) 



with very good convergence, that we do not show again. 
The fields a, f, U and V are shown in Fig. ITH 



B. All other continuations 

We now consider the other possible continuations, in 
particular those of Figs. [21 and [21 Before performing a 
numerical search of the possibilities in the U t space, we 
study the equations of motion in the future patch. We 
proceed in four steps: 

a. A necessary condition for another SSH is a 2 > 2. 
In order to obtain Fig. [3J or any even more exotic con- 
tinuation, we must have a self-similarity horizon before 
the central singularity occurs at x = x s = 0. 

A self-similarity horizon is a DSS line x = Xh (r) (peri- 
odic) where A — 0. The only factor in A = —4a 2 f(f + x) 
that can vanish (while the metric is regular) is / + x. We 
have 



(f + x),x 



(« 2 - 1)/ 



and so 



{f + x), x \ f+x __ 



= 2 



(76) 



(77) 



At the Cauchy horizon / + x — and (f + x) tX ~ 1 
because a 2 — 1 is small, and so (/ + x) > and A < at 
least to the immediate future of the Cauchy horizon. If 
there are more self-similarity horizons to the future of the 
Cauchy horizon, we must have f + x — again there, and 
therefore (/ + x) tX < in some intermediate region, a 2 
must increase from a 2 ~ 1 to a 2 > 2 in order to achieve 
this. 

b. Once a 2 < 1 for any t, a 2 — > for that r. a is 
given by the constraint which reduces at the Cauchy 
horizon to 



(a- 2 ). T = (1 + 2Vf)a 



1, 



(78) 



and this means that 1 < a < 1 + e there. Now a obeys 
the evolution equation 



2x{lna) jX = a 2 - 1 - 2U 2 



(79) 



Recall that in the future patch x < and increases to- 
wards the future. Therefore a — 1 is a repeller in the 
absence of matter, but any outgoing radiation U drives 
a to smaller values. Hence, once a becomes smaller than 
1 for some value of r, it will afterwards decrease to for 
that r whatever happens. Because [i = 2m/r = 1 — a -2 , 
this means that we are in a negative mass regime and 
remain there until we reach the central singularity x = 0, 
which is therefore timelike and has negative infinite mass 
at least for this r. The singularity cannot be reached 
with a finite value of < a 2 < 1. 
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FIG. 13: Variation of the results of the outer patch under changes of the input data from the past patch. On the left we 
represent ||<5/||oo/<5pfc where Sf is the change in 10 s e (continuous line) or in one of the functions W 2 Vf (dotted line), £ (dashed 
line) or U t (dot-dashed line), under a change Spt in the real part of the mode k of U p . On the right we show the same for the 
imaginary part of the modes of U p . We lack enough precision to calculate the last point of the Vf curves. 




-0.9 -0.99 -0.999 -0.9999-0.99999 -0.9 -0.99 -0.999 -0.9999-0.99999 

•^left bright 



FIG. 14: Convergence of the free data on the outer patch. These plots display norms of differences between consecutive 
results: oo-norm (continuous line) and 2-norm (dotted line). On the left, convergence with respect to a;i e ft shows order 3.0 for 
[7 e (triangles) and £ (diamonds) and order 2.4 for Vf (stars). On the right, convergence with respect to alight shows order 2.0 
for £ and Vf and order 1.0 for U t . 



c. Is a 2 > 1 at the singularity possible at least for 
some values of t? If we choose U e much larger than 
0(e 1 ^ 2 ), the U 2 term will drive a to a < 1 almost imme- 
diately. However, we know that the regular solution does 
not have negative mass and therefore it seems plausible 
that functions U e close to the regular case could gener- 
ate a singularity with positive mass. We shall now argue 
that this is not true. 

a can only increase from its CH value if \U e \ is very 
small. In order to explore this regime we return to the ap- 
proximation of perturbing around Minkowski spacetime, 
but now dropping the assumption that F(t) = — G(r). 
The result can be summarized as: 

a(r, x) = 1 + e x 2 a rcg (r, x) 

- ex- 2 [F(r) + G(r)} 2 (80) 
+ eO(x- v ) + 0(e 3 / 2 ). 

The function a rcg is positive and always smaller than 2, 
and is independent of F. All singular terms vanish for 
F = — G. For every other F, the function a becomes 
smaller than 1, which corresponds to negative mass, as 



the center x = is approached. This is due to the di- 
vergent terms (F + G) /x in l|73|) and l|74|) . The crucial 
point is that a 2 is integrated from — U 2 and therefore 
the term in the perturbation expansion that makes the 
mass nonzero has coefficient — (F + G) 2 , and so the mass 
cannot be positive. The negative mass regime is reached 
while perturbation theory still applies, and we have seen 
that afterwards the mass must decrease indefinitely. 

With these arguments we have ruled out the possibility 
of having positive mass at the singularity for very large 
U e and for very small (in particular, close to the regular 
case) U e . However we could be missing some interme- 
diate regime where perturbation theory does not apply. 
For these intermediate cases we observe that the mini- 
mum of a is far from 1 (and therefore we cannot apply 
perturbation theory), but its maximum is positive and 
close to 1, and only becomes negative in the vicinity of 
the singularity. There could be solutions where a is pos- 
itive at the singularity for some values of r. 

Numerically, however, we do not find that. We have 
performed a large number of numerical evolutions in this 
intermediate regime starting from U e — 2" J7 erog (r + 
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FIG. 15: Convergence of V f . Left figure: log 10 |V> S fc-i| for 15 values xieft =-0.8976, -0.9488, -0.9744, -0.9872, -0.9936, -0.9968, 
-0.9984, -0.9992, -0.9996, -0.9998, -0.9999, -0.99995, -0.999975, -0.9999875, corresponding from top to bottom, respectively. It 
seems that all but the first two modes converge to amplitudes below our error threshold when xi e ft approaches -1, even though 
high-frequencies become unstable for the last values of xi e ft. Right figure: modulus of the Fourier transform of consecutive 
differences in the left figure. Convergence is very clear in every mode (including the first one). 
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FIG. 16: Average of V 2 on r for different evolutions from the 
same Ki e ft = —0.9999 to 15 different values of a; r i g ht (those of 
table UTTt . The agreement between the solutions with small 
bright and those with x r ight very close to Xf = 1 implies that 
the final almost constant value is neither a numerical artifact 
nor a consequence of our ansatz. 



toA/8), for n = —5, 10 and m — 0, 7. Even though 
it is difficult to evolve the system near the singularity, 
we always observe a final decay of a to 0. See Fig. ITH1 
Therefore we believe that, starting from small null data 
Vf at the CH it is not possible to form another SSH. 
Therefore, either we have a regular center or a timelike 
singularity, and this singularity always has negative mass. 
In the approach to the center with a 2 — > we always have 
U + V — > with both U, V finite. It should be stressed 
that this scenario is a consequence of the small ampli- 
tude of the null data Vf on the Cauchy horizon in the 
Choptuik solution. We have performed other evolutions 
starting from large null data on the Cauchy horizon (not 
Vf of the Choptuik solution) where a 2 increases beyond 
2. 

d. How is the singularity approached? We find that 
in all numerical continuations a — > 0, / — > oo, U and V 
tend to finite values and U — > — V as the singularity x — 



is approached. If we assume that a — > and U — > Uo(t) 
as x — ► then the a tX and f tX equations to leading order 
give 



/ ~ fo(r)x-\ 



(81) 



a 2 ~ a 2 (r)\x\ 2 ^\ a(r) + C/„(r) 2 . (82) 

Making the ansatz that 

U ~ U q (t) + U 2 (t)x 2 

+ U 21 (t) x 2 In |z| + U a (T) M 2q(t) , (83) 
V ~ V (t) + V 2 (t)x 2 

+ V 2l (T)x 2 ln\x\ + V a (T)\x\ 2a ^, (84) 

we find from the U tX and V x equations to the leading 
order, 0(ln \x\) that 



^0 



-Uq, 



v 2 = U 2 + -fo 1 U^ 



21 



v 2l 



/o Uq, 



-V a = -^-a U . 



(85) 
(86) 

(87) 

(88) 



(The expansion also holds in the special case a = 1.) The 
next order, O(l), gives 



\n(a 2 )' + 1 + 2Ul + 8f U U 2 + 2U Q U' Q = 0. 



(89) 



The metric in the future patch contains a residual gauge 
freedom worth one periodic function of r. Near the sin- 
gularity x — 0, we can fix this gauge freedom by setting 
/o(t) to an arbitrary value. (In our numerical evolutions 
starting at the Cauchy horizon, the gauge has of course 
been fixed already by setting / = 1 at the Cauchy hori- 
zon.) Uq(t) and U 2 (t) are then physical free data which 
determine oo(t) and all other coefficients of the expan- 
sion. (Eq. (|89|l can be solved for an only if the right-hand 
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FIG. 17: Future patch fields on a single A-period for the (probably unique) continuation with a regular timelike center. Note 
that a and / are so close to their flat spacetime values that we plot their difference from 1, rather than the fields themselves. 



side has vanishing average, which means that fo cannot 
be set completely freely.) This expansion is therefore 
generic in the sense of depending on two free functions 
after the gauge has been fixed. 

The behavior just described is what we observe nu- 
merically for large values of the free data U € (t) at the 
Cauchy horizon. For small values of U t {r) we find 
U (t) ~ ±1/a/2, that is a(r) ~ 1. Note that by our 
ansatz of exact DSS with k = 0, U must be an odd func- 
tion of t with zero average. We observe in fact that 
U(t 7 x) goes to a fundamental frequency square wave 
of amplitude 1/V2, that is, in half of each r-period 
U -> 1/V2, and U -> -l/y/2 in the other half. This 
is shown in Fig. ^5] 

As the center x = is approached with Uq = ± 1/ a/2, 
we observe empirically that the r-derivatives become dy- 
namically negligible. This means that different values 
of r effectively decouple, that at each point r the evo- 
lution equations become an ODE system in x, while 
the constraint becomes algebraic, and that the space- 
time becomes locally CSS. It also means that the evolu- 
tion equations become "velocity-dominated" in the sense 
that all derivatives transversal to the singularity (here 
in spherical symmetry, this is only the r-derivative) be- 



come dynamically irrelevant compared to the one deriva- 
tive running into the singularity (here, the rr-derivative) . 
It is known that generic spacelike singularities in general 
relativity with massless scalar field matter are velocity- 
dominated ^(j. Here we find this to be the case only in 
the limit of small data U ( (t) . 

As this class of continuations seems to be locally CSS 
near the singularity, it is interesting to study the exactly 
CSS solutions from the point of view of a DSS ansatz. 
Starting from a generic DSS scalar field (|12fl we introduce 
a new (gauge-dependent) variable 



W = U + —(U + V) = V2nG(n + ib T ), (90) 

X 

which coincides with — V at the future lightcone j jx — 
— 1, and obeys the equation 

xW, a = -U T . (91) 

In exactly CSS solutions of the system the metric func- 
tions and U, V are independent of r, and hence W is con- 
stant with value W = \J2-kGk. The CSS solution with 
k = was found in closed form by Roberts and is 
described in Appendix [U] The CSS solutions with k ^O 
were studied numerically by Brady |12| . 
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FIG. 18: Maximum and minimum of a for different initial 
data functions U t from numerical evolutions of the nonlinear 
equations. Note that the x axis becomes logarithmic when 
approaching the boundaries, while the a axis is linear between 
the dotted lines at a = 1 — 10~ 6 , a = 1 + 10~ 6 , and logarithmic 
elsewhere. The initial U ( have been obtained multiplying the 
regular data (black) by different constants: (red), -1 (blue) 
and 5 (green). Clearly, the larger [7 e is, the sooner a decays. 

We find that the small t/ e (r) continuations locally ap- 
proach a Roberts solution (with n = 0) with a value of 
the parameter p of the Roberts solution that depends on 
r, and not, as one might have expected, one of Brady's 
solutions. This is discussed in Appendix IU1 

VI. GLOBAL IMAGES OF THE CHOPTUIK 
SPACETIME 

Figures 03 ^1 \^\ an d El show in a very detailed way 
the structure of the Choptuik spacetime. However it is 
difficult to get from them an idea of what it looks like 
globally. In this final Section before the conclusions we 
present a number of additional figures that will fill this 
gap. We do it in two very different ways. 

A. Global coordinate systems 

As shown in Fig. ^ our three (r, x) patches match 
continuously, but we do not expect the resulting global 
coordinate system to be differentiable at the interfaces 
between the patches. The critical spacetime itself, how- 
ever, is differentiable (analytic at the past lightcone and 
C 1 at the CH) , and it must be possible to construct global 
coordinate systems which are at least C 2 . 

One simple possibility is synchronous slicing plus area 
gauge: 

ds 2 = -dT 2 + 2B(T, R)dTdR + C(T, R)dR 2 + R 2 dSl 2 . 

(92) 

We add the gauge condition T = —R on the past light- 
cone of the singularity. The coordinate transformation 
T(t,x), R(t,x) can be easily integrated and it is shown 
on the left panel of Fig. [20|for a single A-period in r. 



Another simple possibility is double null coordinates: 

ds 2 = -oj(u, v)dudv + R 2 (u, v)dn 2 (93) 

with gauge condition u = v = T at the center, where T 
is the time coordinate constructed in the previous para- 
graph. The coordinate change u(t, x), v(t, x) is shown on 
the right panel of Fig. [20|for a single A-period in r. 

The fields a, U and V in double null coordinates are 
shown in Fig. [2] As these fields are spacetime scalars, 
they should be analytic in these coordinates at the past 
lightcone, and C 1 at the Cauchy horizon. In the plot it 
looks as if they are only C° at the Cauchy horizon, but 
that is due to a lack of resolution in the plot: the slopes 
on the two sides of the Cauchy horizon are dominated by 
ai+ e and ai+2e, which are discontinuous at the resolution 
of the plot, but close enough to the Cauchy horizon, the 
slope becomes ai, which is continuous. 

B. Dynamical phase space portraits 

We shall now consider the spherical DSS scalar field 
as an infinite-dimensional dynamical system where x is 
the "time" coordinate. The dynamical variables in this 
system are U(t), V(t) and f(r) (or b(r) in the outer 
patch). The variable a(r) is not independent, but given 
by a constraint. However, many solutions of the dynami- 
cal system correspond to the same spacetime, namely all 
that are related by the coordinate transformations JSJ. 
The pair of periodic functions U(t) and V(t) describes 
a closed, possibly self-intersecting curve in the (U, V) 
plane. An entire evolution in x gives a surface that is 
topologically a cylinder, which we may call a phase por- 
trait of the solution. Clearly the surface itself is invari- 
ant under the coordinate change r — > r + ip(x, r). If we 
consider r as our "space" coordinate in the usual "3+1" 
split, then by looking at the phase portrait we have elim- 
inated the "spatial" gauge freedom. The "slicing" free- 
dom x — ► tp(x, t) however does change the shape of the 
phase portrait, so it is not completely gauge-invariant. 
The Choptuik spacetime up to the CH is given in Fig. 1221 
and its regular continuation is shown in Fig. 1231 

Imposing CSS means that the system is independent 
of t and hence the phase-portrait in (U, V, x) reduces to 
a line that can be easily projected on the (U, V) plane. 
Then the whole evolution of the system can be described 
by this curve in the (U, V) plane, which is now completely 
gauge- invariant . This is essentially what Brady has done 

In order to understand the singular continuations of 
the critical spacetime, we first look at phase portraits 
of the Roberts spacetime. Fig. 1241 shows the phase-lines 
of the Roberts solution for several values of p in both 
branches (see also Appendix 0. We see that the shape 
of the curves become constant for very small values of 
p, but just translated along the logx-axis. On the other 
hand Fig. 1251 gives the phase-portrait two of our singu- 
lar evolutions, for small and large initial U e , respectively. 
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FIG. 19: Future patch fields on a single A-period for a singular continuation of the Choptuik spacetime. 



In both cases as U + V — > squeezes the phase cylin- 
der into the diagonal, following a Roberts solution in the 
former case, and not doing so in the latter. This figure 
demonstrates why it is not possible to obtain a 2 = 2 in 
the continuations of the Choptuik solution. The Roberts 
spacetime with p = 1 does have a 2 — 2 on the singular- 
ity, p = (Minkowski) has a 2 = 1, and all other p have 
a 2 — > 0. The p = 1 Roberts solution has the phase line 
U = V = ±l/y/2 line. To reach a 2 = 2 approach a so- 
lution must approach this line from the very beginning, 
because it is unstable. This requires V of order 1/V2 at 
least and on the Cauchy horizon we only have V of order 
10" 3 . 



VII. CONCLUSIONS 

In a companion paper Q we investigate the con- 
straints that the kinematic assumptions of spherical sym- 
metry and discrete self-similarity impose on the causal 
structure. The key elements of our analysis were the 
self-similarity horizons - radial null geodesies that are 
mapped onto themselves by the self-similarity. These 
come in two types: a "fan" connects a point singular- 



ity to a piece of null infinity, while a "splash" connects 
a piece of null singularity to a point at null infinity. All 
spherically symmetric and self-similar spacetimes can be 
enumerated in terms of a sequence of fans and splashes. 
We found that the singularity is central and consists of 
a middle segment which is either a point or null line, 
flanked by two segments which can be timelike, space- 
like, or null lines, or which can be absent. There are 
never two or more disconnected singularities. 

In this paper, we have focussed on scalar field matter 
and in particular Choptuik's critical solution in gravita- 
tional collapse. This solution is found as an intermediate 
attractor in the evolution of regular asymptotically flat 
scalar field initial data near the black hole threshold pj. 
Independently, it can be constructed as the solution of 
a nonlinear PDE boundary value problem from the as- 
sumptions of spherical symmetry, discrete self-similarity, 
and analyticity at the past center and past lightcone of 
the singularity 5]. It can then be uniquely continued up 
to the future lightcone (Cauchy horizon) of the singular- 
ity. Ref. [| found that the scalar field oscillates roughly 
as cos (In \y\) as the Cauchy horizon y = is approached. 
It was suggested on theoretical grounds that these oscil- 
lations are damped by a factor \y\ £ with e > 0, but the 
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FIG. 20: Spacetime diagrams for a single period of the Choptuik spacetime in differentiable coordinate systems (synchronous- 
area on the left and double null on the right). It is possible to construct the whole spacetime by adding rescaled (by a factor 
e A = 31.357 or its inverse) blocks to the center and outside the figures. Continuous lines represent r = const, lines and dashed 
lines represent x = const, lines. By our gauge choices times T and (v + lt)/2 coincide on the centre worldline. 



numerical value of e was too small to distinguish it from 
zero numerically. 

We have repeated the numerical calculations of Ref. p| 
from scratch and have increased the accuracy by more 
than four orders of magnitude. We have found the pos- 
itive value for e given in Eq. (|63(l . Therefore the scalar 
field is continuous on the Cauchy horizon, with null data 
given in Eq. (J7DJ. The structure of the fields near the 
Cauchy horizon is quite complicated, and has been dis- 
cussed in Subsections IIII Cl and IIII Dl 

All possible DSS continuations beyond the Cauchy 
horizon are determined by the (fixed) null data on the 
Cauchy horizon, and one (arbitrary) periodic function 
U e (t) that can be thought of as data emerging from the 
naked singularity. (In the absence of the naked singu- 
larity, null data on the Cauchy horizon would alone de- 
termine the continuation.) There is a unique DSS con- 
tinuation that has a regular timelike center except for a 



naked point singularity. For all other values of the ar- 
bitrary function U e (f), the continuation has a timelike 
central singularity with infinite negative mass. In partic- 
ular, neither Fig. |3 nor the even more exotic spacetime 
diagrams found in self-similar perfect fluid solutions Q 
can arise as continuations of the Choptuik solution (al- 
though they can probably arise in other DSS scalar field 
solutions). 

The global structure of the Choptuik solution is of in- 
terest partly because of the connection between critical 
collapse and cosmic censorship. Choptuik's work has es- 
tablished that (assuming the validity of general relativity 
at arbitrarily high curvatures) a naked singularity can be 
formed in the spherical Einstein-scalar field system from 
generic regular and asymptotically flat initial data by 
fine-tuning any one parameter to the black hole thresh- 
old. Solutions with a naked singularity therefore form a 
subset of codimension one (the black hole threshold in 
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FIG. 21: Scalars a, U and V as functions of the double null variables u and v. Grid lines are lines of constant x or r. Time is 
increasing from left to right, and the central world line is at the back. 



the space of initial data) in all solutions arising from reg- 
ular data. The Choptuik spacetime is an attractor on 
the critical surface and therefore, assuming that it is ac- 
tually a global attractor on the surface, we can conclude 
that every naked singularity will have, at least in a neigh- 
borhood of the singularity, the structure of the Choptuik 
spacetime. 

We can now summarize this structure as follows: The 
curvature at the Cauchy horizon of the singularity is fi- 
nite but not differentiable, with an infinite number of 
damped oscillations piling up as the Cauchy horizon is 
approached. The continuation of the spacetime beyond 
the Cauchy horizon is not unique, but we have shown 
that if it is DSS then the naked singularity is either a 
single point or timelike with infinite negative mass. The 
spacetime near the timelike singularity is locally CSS and 
velocity-dominated . 

The question remains whether this structure is stable 
against perturbations which break self-similarity and/or 



spherical symmetry. Because the background is spheri- 
cally symmetric and periodic in t, all such perturbations 
are of the form e Ar f(x, r), where f(x, r) is periodic, times 
a suitable scalar, vector or tensor spherical harmonic. We 
have shown in the past 0, 0] that all but one of these 
modes decay, with ReA < 0. (The one growing mode de- 
termines the critical exponent for the black hole mass.) 
The functions f(x, r) were calculated explicitly only be- 
tween the regular center and past lightcone, but as they 
are by construction analytic at the past lightcone, they 
can be analytically continued up to the Cauchy horizon, 
with the same A. 

What needs to be done is to check that the functions 
f(x, t) remain bounded as x approaches the Cauchy hori- 
zon. We expect that this is true. Nolan and Waters 
|14| have investigated a massless minimally coupled, non- 
spherically symmetric scalar field propagating on a class 
of spherically symmetric CSS spacetimes, and find that 
its gradient remains finite at the Cauchy horizon. The 
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FIG. 22: Left: Phase-portrait of the Choptuik spacetime from the regular past centre at x c to the CH at Xf. The a;-axis has 
been highly distorted in order to show all interesting details: the axis is logarithmic between x c and x p with accumulation point 
at x c (note that the label x c has been situated at finite distance for convenience); the axis is logarithmic with accumulation at 
Xf between x v and xf, the axis has been transformed from x to Xf — x 1 '* to show the decay of the function U towards Xf (this 
is a semi-analytical extrapolation of our numerical data). Vertical generatrices are r = const, lines. Right: a reduction in the 
V axis in the region between x v and Xf to show the oscillations. Compare with Fig. 1121 




FIG. 23: Phase-portrait of the regular continuation of the 
Choptuik spacetime from the CH at Xf = —1 to the regular 
center at x T = 0. Again, from x — — 1 to x = —0.999 we have 
distorted the axis to show the infinite number of oscillations 
that pile up there. From x = —0.999 to x = —0.5 we use a 
logarithmic axis and from this middle point to x = —0.001 
we use a different logarithmic axis. 



structure of the full perturbation equations is similar. 

As the functions f(x, t) cannot be analytic at the 
Cauchy horizon (because the background is not), the 
spectrum of A need not be the same in the continua- 
tion beyond the Cauchy horizon, and is likely to depend 



on the choice of continuation. 

We have finally given, for the first time, a number of 
global images of the Choptuik spacetime, trying to min- 
imize the gauge content of the pictures. 
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APPENDIX A: REVIEW OF COORDINATE 
SYSTEMS FOR GIVEN F 

In this Appendix we review a variety of coordinate 
systems defined by giving F as a function of r and x. 
This class covers the coordinate systems used in Refs. 
[E S 113 El j but does not include the Lagrangian 
fluid coordinates used in ^tJ- We classify the gauges 
within this class in three stages. 1) We fix F(t,x). 2) 
We impose an algebraic relation between A, B and C . 3) 
We parameterize A, B and C in terms of two functions, 
say a and /. (It will turn out to be useful to always use 
the scalar a defined above as one of the two parameters.) 

In a DSS spacetime, the unknowns are periodic in r. 
We therefore solve the Einstein equations by "evolving" 
in x, with periodic boundary conditions in r. When we 
use the Einstein equations to eliminate the metric deriva- 
tive P, U and V obey a pair of transport equation of the 
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FIG. 24: Phase-lines of the Roberts spacetime. The colors 
encode different values of p, for both branches: 1 (red), 0.99 
(orange), 0.9 (yellow), 0.5 (light green), 0.1 (dark green), 0.01 
(light blue), 0.001 (dark blue) and 0.0001 (purple). The future 
(past) branch of the singularity is denoted with thick (thin) 
lines starting with V = (U = 0). Note that the vertical lines 
U — V — ±l/\/2 corresponding to p = 1 are unstable. Note 
also that the smaller p is, the longer the curves stay close to 
the unstable flat line U = V = 0. 

form 

U, x = {...)U, r + (...), V x = (...)V r + (...), (Al) 

where the dots stand for a known function of U , V, a and 
/ but not their derivatives. By a suitable parameteriza- 
tion of A, B and C in terms of a and /, the Einstein 
equations can always be brought into the form of two 
ODEs in x, which for our purposes are evolution equa- 
tions, and one ODE in r, which for our purposes is a 
constraint: 

/,»=(...). o,x = (•••), a,r = (•••)• (A2) 

This last equation can be made linear in a~ 2 (that is, in 
/j,) by a suitable choice of /, at least in all three coordinate 
systems we use. This inhomogeneous linear equation can 
then be solved uniquely for a in terms of U, V and /, 
using the fact that U , V and / are periodic in r and 
that we require a to be periodic, too. In a CSS solution, 
where nothing depends on r, the a_ T constraint becomes 
an algebraic equation linking U , V, a and /. 

1. Past patches with F = x 

The regular center is both timelike and an x line, but 
r is finite there. This requires F = at the center, and 
F > elsewhere. The simplest choice is F = x. On 




FIG. 25: Phase-portrait of two singular continuation of the 
Choptuik spacetime (upper: initial U e of order 10 -3 ; lower: 
initial U e of order 10), together with a few Roberts solu- 
tions for p = 0.1,0.01,0.001. The p = 0.001 solution al- 
most coincide with the external envelope of the upper phase- 
portrait and sets a maximum for the values that are realized 
in this solution. The lower phase-portrait is shown only for 
—0.999 < x < —0.4. For —0.4 < x < it is essentially on the 
diagonal in the (U, V) plane although the amplitude of the 
oscillations is not the same as in the Roberts solution. Note 
that for small x this phase portrait has been truncated also in 
U, because U becomes large. For x ~ — 1 the phase portrait 
is essentially on the ?7-axis. 

a past patch, it is possible to choose the r lines to be 
timelike, null or spacelike, or to change signature. We 
consider only the first two possibilities 

e. Bondi past patch If we choose the r-lines to be 
null everywhere, this means C = 0. We then have B < 
between the past center and the past light cone, and 
B > beyond the past lightcone, so that A and B both 
change sign at the past lightcone. F = x and C = has 
been used by Brady 12] to investigate CSS solutions with 
scalar field matter. He used the Bondi metric coefficients 
that are traditionally called g and g as parameters. In 
a DSS solution this gives a constraint for a combination 
of g and g, as well as evolution equations for g and g. 
The best parameterization we have found uses the mass 
function a and the Bondi metric coefficient / = g/a 2 : 



A = -g(gT2x) = -a 2 f(fT2x), (A3) 

B = T g = T a 2 f, (A4) 

C = 0, (A5) 

F = x, (A6) 



where the upper sign applies when r is an outgoing null 
coordinate, and the lower sign applies when r is an ingo- 
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ing null coordinate, assuming that / > 0. This gives an 
ODE evolution equation and a linear ODE constraint for 
a -2 , and an ODE evolution equation for /. 

/. Schwarzs child past patch If the r-lines are space- 
like everywhere, this implies C > 1. A natural algebraic 
condition to impose is to make the r lines orthogonal to 
the r lines. This means B(r,x) = —xC{t,x). t is then 
— ln(— t) where t is the Schwarzschild time coordinate. A 
possible parameterization is in terms of the metric coeffi- 
cients often called a and a in Schwarzschild coordinates, 
(a is the scalar we have defined above) . The best param- 
eterization we have found is by a and / = a/ a, and is 
given in above. It gives ODE evolution equations for 
a and /, and a linear ODE constraint equation for a~ 2 . 
The transport equations for U and V are linear in U and 
V because P reduces to a function of a, t and x. The 
lightcones are at / ± x = (although of course only one 
lightcone can be covered at one time by a past patch) . 



2. Outer patches with F = 1 

On an outer patch that stretches from the past to the 
future lightcone the r lines must be timelike at least 
somewhere, and it is possible to make them timelike ev- 
erywhere. The simplest choice for F compatible with all 
these possibilities is F=l. 

g. Schwarzschild outer patch Assuming that the r 
lines are timelike everywhere means C < 0. We have B < 
on the past lightcone and B > on the future lightcone, 
so that B has to change sign somewhere between. This 
"B-surface" is where the r and x lines are orthogonal. 
In a one of us used a coordinate system in the outer 
patch that was based on Schwarzschild coordinates. As in 
the past patch based on Schwarzschild coordinates, this 
meant B = —xC. However, when we impose B = — xC 
together with F = 1, we face an unexpected coordinate 
singularity at x = 0. In particular, if we again use a and 
f = a/a to parameterize the metric, we obtain an ODE 
constraint and an ODE evolution equation for a as before 
on the past patch, but for / we obtain an equation of the 
form 



f,x 



(...)/,r + (■■•) 



(A7) 



In order to make the coordinates regular at x = 0, we 
need to impose the vanishing of the numerator of this 
equation there, and that is effectively what was done in 
[5j . But this introduces an additional boundary condition 
just to keep the coordinate system regular on a surface 
where the solution itself is perfectly regular, and that is 
why we do not use these coordinates here. 

h. Best buy outer patch A better algebraic condi- 
tion to combine with F = 1 is C = — 1. A workable 
parameterization is in terms of the mass function a and 
B itself. This gives ODE evolution equations for a and 
B, and a (nonlinear) ODE constraint for a. The light- 
cones are at a ± B = 0. If we replace B by b — B/a, we 



have the parameterization given above in i|21[l . and the 
ODE constraint for a becomes linear in a~ 2 , while the 
wave equation becomes linear in U and V . 

i. A — —C: unworkable It is compatible with F = 1 
to assume that the t lines "bend round" to become space- 
like at large and small x. means This that C changes sign. 
A natural choice is that A and C change sign together, 
and we can impose this by setting A — — C. We have 
not found a way of parameterizing this choice in a way 
that brings the Einstein equations into the standard form 
(|A2|) . Parameterizing with B and fi gives two roots for 
A. Using A and /i gives two roots for B. 

j. x and t orthogonal: unworkable We can also force 
the "bending round" by making r and x lines orthogo- 
nal everywhere, or B — 0. This implies that, while they 
are independent, A and C change sign together. Param- 
eterizing by A and C we obtain an evolution equation 
and a constraint for A and a constraint for C. But the 
constraint for C is homogeneous, C )T = N(A, U, V)C, 
so that C cannot be periodic. Therefore this coordinate 
system is not compatible with CSS or DSS. 

k. Forcing the lightcones: unworkable Instead of 
finding a working outer patch and then forcing the two 
lightcones to fall on the lines x = x p and x = Xf by im- 
posing A(t, x p ) — A(r,Xf) — 0, we can make A a given 
function of x. The simplest such choice is A — 1 — x 2 , 
which vanishes at x = ±1. (The two lightcones will be 
distinguished by the sign of B). A parameterization that 
brings the Einstein equations into standard form is 



,4=l-a; 2 , B = bc, C 



F = 1. 



(A8) 



This gives an evolution equation and a constraint for b 
and a linear constraint for c. However, the right-hand 
sides of all three equations are of the form N/b. Fur- 
thermore, at b = the numerator of 6 r depends on the 
matter through U 2 + V 2 , while the numerator of fe x de- 
pends on the matter through U 2 — V 2 . (The numerator 
of c r is proportional to that of b , T .) That means that we 
would have to impose two independent regularity con- 
ditions at b — 0. Together these would fix U and V 
completely. Therefore this coordinate system is not suf- 
ficiently generic near the line 5 = 0. 



APPENDIX B: SINGULAR EXPANSION 
AROUND THE CAUCHY HORIZON 



From I|46I47|) we see that 

D = yD 1 (r)+y 2 D 2 (T)+0(\y\ 2 ^) 



(Bl) 



The regular coefficients to 0(y) are as follows: U\(t) is 
the unique periodic solution of 



Ui + (1 - al - D x )Ui + Ui - 2a oai U = 0, (B2) 



and 



Vi - -^-[(l-a 2 )V + U + Vj], 



(B3) 
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a i 



-ku 2 -v % 



h = -£--^(-3 + a 2 + 2V Q 2 ), 
Zao 



(B4) 
(B5) 



-J- (2a oai - b ± U 2 + hV Q 2 + WoVx) , (B6) 
aoh (B?) 



The regular coefficients to higher orders continue in this 
style, with U n {r) the solution of a linear inhomogeneous 
ODE, while V n , a n and b n are given algebraically. 

The nonvanishing singular coefficients up to 0(y 1+3e ) 
are 

a 1+e (r,x) = a 1+e (r)ai +£ (f), (B8) 

ai +2 e(T,x) = ai +2 e{T)a 1+2e (T), (B9) 

V 1+€ (t,x) = V 1+£ (t)V 1+€ (t), (BIO) 
U 1+e (r,x) = TMU-rWiUr), (Bll) 



i=i 



[/ 1+2e (r,x) = ^^ 2e (r)^ 2e (f), (B12) 



f/ 1+3e (r,x) = ^U^ 3e (r)U[% e (r). (B13) 



We only give two examples of how these coefficients are 
derived. Substituting the ansatz into the a tX equation 
and isolating the terms of 0(y e ) in the result, we obtain 

(1 + e)ai +£ (r)ai +e (f) + Ka 1+e (r)a' 1+e (f) = 

-ar)U (T)U e (T)U e (r). (B14) 

We solve this for all r and f by setting 



ai+eO) 



gr)U (T)U e (T) 
K 



(B15) 



and by making di +e (f) the unique solution of the ODE 
1 + e. 



K 



-a 1+e + U e = 0. 



(B16) 



If we substitute the ansatz into the U x equation l|27|) 
and isolate the terms of 0(y 1+e ), we find 



= (1 - «o) Y, U ilMle ~ 2a oai U c U e - 2a a 1+e a 1+e U + V 1+e V 1+e + ^ U'^Ml* + (1 + E ^lMS ( B17 ) 

Because KD\ = 1 + H', the derivatives E^+ e cancel out. Taking also into account that Vi+ e ~ ai +e (this is an accident 
at this particular order), the equation can be rewritten as 



£ + (1 - og - (1 + ejA) U[ 



(0 



[-cD a - 2a a 2 )tfe] £> £ + [-iO? 2 fr e ] ^ + [Vi +e - 2a ai+eC/o] ai+ e - 0. (B18) 

I 



The terms in square brackets all depend on r and each 
multiplies a different known function of f . To solve this 

* (i) 

for all r and f , we assign one term U 1 f e to each function 
of f : 

U[% = U € , U[% = U' e , U[% = ax +e . (B19) 

The corresponding coefficients u[ l ^_ e are the unique solu- 
tions of the ODEs 

U&l + (1 - al - (1 + e)D x ) tj[% + S« = 0, (B20) 

where the source terms can be read off directly from 
(|B18f) . The calculation of the other terms proceeds in a 



similar manner at all orders. Note that the functions of 
type /, where / stands for V, a and 6, are given alge- 
braically and / obeys an ODE. For U it is the other way 
around. 

The limit Vq = e = of our expansion exists. 
The regular 0(y n ) terms in the series vanish identically, 
and so do some of the 0(\y\ n+ke terms but not all of 
them. eti and eti +e for example vanish because they are 
proportional to Uq, but a x+ 2e is proportional to U 2 , and 
so encodes the curvature to 0(y). 

In the limit e = the curvature components propor- 
tional to U 2 are no longer continuous at the CH because 
they are now periodic in f = r — In \y\, but for the same 
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reason they remain bounded. The components propor- 
tional to UV and V 2 are still continuous even then be- 
cause V is 0(y). a is also still continuous, with a = 1 on 
the Cauchy horizon. 



APPENDIX C: THE ROBERTS SOLUTION 

In the notation of [lg, the Roberts solution is 
given by 



ds 2 



= —dudv + r 2 (u,v)dQ 2 , 



r 2 (u,v) — — [(1 — p 2 )v 2 — 2vu + u 2 ~\ , 



(f>(u,v) = 



log 



(1 — p)v — u 



(CI) 
(C2) 

(C3) 



with p a constant parameter, p — is Minkowski space- 
time with zero scalar field, and without loss of generality 
p > 0. Only the regions r 2 > are physical and without 
loss of generality, we consider the right side (v — u > 0) 
of the spacetime. 

For p ^ 0, the lines u = (1 ± p)v are central curva- 
ture singularities: the mass function /i = —p 2 uv/(4r 2 ) 
diverges on them. The line u = (1 + p)v, v < Ois 
timelike and has (infinite) negative mass, and for v < 
forms the past branch of the Roberts singularity. The 
line u = (l—p)v is timelike with negative mass for p < 1, 
null with zero mass for p = 1, and spacelike with positive 
mass p > 1. For v > it forms the future segment of the 
singularity. 

Some of our solutions become asymptotically CSS and 
have k ~ 0, therefore approaching the Roberts space- 
time. We know that the solutions develop a timelike, 
negative mass singularity, which must correspond to one 
of the two branches of the Roberts singularity. In this Ap- 
pendix we study the issue of which of the two branches is 
actually approached in our numerical evolutions and for 
what values of p. 

Our scalars U and V are 



U 



P 



V2 v 



V 



P 



y/2 V — U — p 2 V ' 



(C4) 



where both denominators are positive in the region of 
interest. It is clear that p is covariantly defined by U 
and V on the light cones u = or v — 0, respectively. 
However, at the singularities u = (1 ± p)v the scalars 
take values U = —V = ±l/\/2 for any =/= p =^ 1 in the 
Roberts solution. If we use an expansion using a generic 
self-similar coordinate y such that u/v = F(y) with the 
singularity at y — [that is, F(0) — 1 ±p]: 



V2U(y) 
y/2V{y) 



±1 



-F'(0) 
P 



0{y 2 ), 



Tl + yl±P F ' {0) + O(y 2 ). 
p 1 ± p 



(C5) 
(C6) 



Therefore p is covariantly given by the ratio of the rates at 
which U and V approach their values at the singularity: 



,. dU l±p 
hm = . 

x^o dV 1 =F P 



(C7) 



This determines both p and which branch of the singu- 
larity we approach locally: The upper sign applies for 
\dU/dV\ > 1 (past branch, turned upside down) and the 
lower sign for \dU/dV\ < 1 (future branch of the singu- 
larity) . 

We can relate the coordinates (u, v) of (|C1IC3|I to our 
coordinate x through 



(l-^) 2 =P 2 + {1-P 2 )x 2 



(C8) 



This expression applies only for negative x and so only 
covers the future branch of the singularity. We obtain 
the expressions for the past branch exchanging the role 
of the functions U and V. That explains why exchanging 
U and V in Eq. (|C7fl amounts to a change of branch. It 
also shows that we have to change branch four times per 
A-period as we move through the four quadrants deter- 
mined by the lines U + V = and U — V = 0. 
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